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""The  objective  of  this  study  is  to  develop  the  basic  equation  for  a  finite  element  formulation 
which  can  be  used  for  solving  problems  related  to  coupled  thermomechanical  behavior. 

The  formulation  is  based  on  the  introduction  of  a  new  quantity,  defined  as  heat  displacement 
which  is  related  to  temperature  in  the  same  manner  as  the  mechanical  displacement  is  related  to 
strain.  The  introduction  of  such  a  quantity  into  the  thermal  equation  leads  to  a  form  similar  to 
the  equation  of  motion,  so  that,  the  governing  equation  for  heat  conduction  can  be  regarded  as 
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20.  Abstract  (Continued) 

Using  this  definition  for  the  heat  displacement  the  coupled  thermal  and  momentum  equations 
are  written  in  a  unified  form  with  the  thermomechanical  coupling  terms  appearing  explicitly  in  the 
constitutive  relation  and  not  in  the  equation  of  motion  or  heat  conduction.  This  presentation  is  more 
appropriate  since  the  behavior  of  a  material  is  expressed  through  its  constitutive  relation. 

The  governing  equations  are  used  to  write  a  variational  formulation  which,  together  with  the 
concept  of  generalized  coordinates,  yields  a  set  of  differential  equations  with  the  time  as  the  independ 
ent  variable.  These  equations  are  presented  in  a  form  equivalent  to  that  of  the  Lagrangian  equations 
in  mechanics. 

Following  this  formulation,  the  basic  finite  element  approach  is  used  to  derive  two  element 
models,  one  which  is  represented  by  a  linear  approximation  of  the  displacements,  and  a  second  one 
which  is  represented  by  a  higher  order  approximation  of  the  displacements. 

An  application  of  the  derived  formulation  for  the  two  element  models  is  given  for  the  problem 
of  induced  thermal  waves  in  a  medium  which  is  subjected  to  temperature  changes  at  its  boundary 
surface.  Results  obtained  for  this  problem  are  for  the  coupled  and  uncoupled  thermomechanical 
behavior. 

The  results  obtained  by  the  present  formulation  are  compared  with  existing  results  of  analytical 
solutions. 
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FINITE  ELEMENT  MODELING  OF  THE 
COUPLED  THERMAL  AND  MOMENTUM 
DIFFUSION  EQUATIONS 


INTRODUCTION 

The  behavior  of  solids  under  temperature  changes  has  long  been  studied  and  a  large  number  of 
publications  exist  in  both  the  fields  of  thermal  and  mechanical  deformations.  Developments  in  both  of 
these  fields  during  the  last  two  decades  has  produced  a  variety  of  solutions  and  techniques  of  direct 
engineering  application. 

The  basic  theory  for  thermal  effects  on  mechanical  deformations  has  been  well-established  and 
analytical  methods  of  solution  exist  for  the  governing  equations.  However,  the  practical  application  of 
analytical  methods  presents  difficulties  for  bodies  of  complex  geometries  or  bodies  under  complicated 
boundary  conditions  and  cannot  be  relied  upon  for  the  solution  of  any  problems  of  practical  importance. 
Consequently,  numerical  analyses  have  been  extensively  studied  and  various  approaches  based  on  varia¬ 
tional  methods  have  been  developed.  Such  methods  and  solutions  can  be  found  for  a  variety  of  appli¬ 
cations.  A  variational  principle  for  the  equations  of  coupled  thermoelasticity,  introduced  by  Biot,  and 
applications  to  this  variational  formulation  are  given  in  (1),  [2]  and  [3]  Boley,  in  [4]-[6],  refers  to  solu¬ 
tions  of  thermoelastic  problems  and  a  problem  on  thermoelastic  wave  propagation  is  solved  by  Dunn  in 
(19].  The  thermal  stresses  induced  in  a  semi-infinite  body  due  to  suddenly  applied  step  heat  input  are 
discussed  by  Grimado  in  [14].  In  addition,  a  variational  formulation  is  given  by  Nicked  and  Sackman 
in  [IS]  for  the  coupled  linear  thermoelasticity  and  an  application  of  this  formulaton  to  a  boundary  value 
problem  can  be  found  in  [16]. 

The  majority  of  the  numerical  analyses  employ  the  finite  element  method  because  of  its  superior¬ 
ity  over  finite-difference  techniques.  The  finite  element  method  has  been  well-established  in  recent 
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years  and  the  basic  formulation  lor  this  method  is  given  under  the  name  of  stiffness  and  deflection 
analysis.  An  introduction  to  the  finite  element  method  can  be  found  in  the  texts  by  Desai  and  Abel 
[111  and  by  Zienkiewicz  [211,  with  applications  to  engineering  related  subjects.  The  method  has  also 
been  employed  in  solving  thermomechanical  problems  and  applications  of  the  method  are  given  in  [17] 
and  [18].  In  most  of  these  studies,  the  effect  of  thermomechanical  coupling  is  included  in  the  formula¬ 
tion  with  applications  to  various  problems  of  practical  interest. 

Existing  finite  element  solutions  are  based  on  formulations  where  the  coupled  equations  are 
treated  as  two  separate  equations  with  the  coupling  appearing  explicitly  in  the  thermal  diffusion  equa¬ 
tion.  As  a  result,  existing  solutions  usually  involve  a  time-consuming  iterative  process.  A  different 
approach  which  uses  finite -difference  methods  to  treat  the  heat  equation  and  finite  elements  for  the 
momentum  equation  was  also  suggested.  Again,  it  requires  a  lengthy  numerical  iteration.  Obviously,  a 
more  appropriate  formulation  will  be  one  which  describes,  by  a  single  discrete  element,  both  thermal 
and  mechanical  deformations.  This  motivates  the  research  reported  in  this  study.  Essentially,  the  pur¬ 
pose  of  this  study  is  to  formulate  both  the  problems  of  heat  transfer  and  thermal  deformation  under  a 
unified  formulation.  In  order  to  achieve  this,  a  unified  variational  formulation  is  introduced  for  the 
coupled  thermal  momentum  diffusion  equations.  The  key  to  this  is  the  introduction  of  a  new  quantity 
defined  as  heat  displacement,  which  is  equivalent  to  the  mechanical  displacement  and  has  the  units  of 
length.  As  a  result  of  this  definition,  the  change  in  temperature  is  treated  as  a  thermal  deformation  and 
it  is  equivalent  to  a  mechanical  deformation  [22], 

In  the  first  chapter,  the  basic  definitions  are  introduced  and  the  coupled  equations  are  written  in  a 
form  such  that  the  equation  of  motion  and  the  equation  of  equilibrium  for  heat  transfer  are  of  similar 
form,  with  the  coupling  term  appearing  explicitly  in  the  constitutive  relations.  Such  formulations  are 
more  appropriate  since  the  behavior  of  a  material  is  expressed  through  its  constitutive  relations  and  not 
through  the  equation  of  motion  or  equilibrium.  A  variational  equation  is  derived,  based  on  the  princi¬ 
ple  of  virtual  work  in  mechanics,  which,  with  the  use  of  generalized  coordinates,  leads  to  a  unified 


2 


NRL  MEMORANDUM  REPORT  4310 


equation  for  the  general  coupled  case.  Since  this  equation  is  expressed  in  terms  of  generalized  coordi¬ 
nates,  which  can  represent  both  mechanical  and  heat  displacements,  it  can  also  be  used  for  deriving  a 
finite  element  formulation.  This  is  done  in  the  second  chapter,  where  the  basic  finite  element  method 
is  used  to  derive  two  finite  element  models  for  solving  initial/boundary  value  problems.  The  first 
model  is  based  on  a  linear  approximation  of  the  displacements  and  the  second  on  a  higher  order 
approximation.  The  matrix  equation  for  the  first  model  is  expressed  in  terms  of  nodal  displacements 
and  for  the  second  model,  in  terms  of  nodal  displacements  and  nodal  deformation. 

Each  element  model  can  be  used: 

a)  for  coupled  diffusion  problems, 

b)  for  uncoupled  ones  by  setting  the  thermomechanical  coupling  parameter  equal  to  zero, 

c)  for  pure  thermal  diffusion  problems  by  setting  the  mechanical  displacements  equal  to  zero,  and 

d)  for  either  dynamic  or  quasi-static  solutions  by  retaining  or  neglecting  the  inertia  term,  respec¬ 
tively. 

Numerical  solutions  of  the  matrix  differential  equations  are  obtained  through  a  third-order  backward 
finite  difference  scheme.  This  implicit  integration  technique  has  been  found  to  produce  stable  results 
and  its  accuracy  has  been  tested  on  a  variety  of  problems. 

The  derived  finite  element  formulation  is  used  in  the  last  chapter  to  solve  an  inital/ boundary 
value  problem  of  transient  thermomechanical  behavior,  when  uniform  heating  is  applied  on  the  boun¬ 
dary  of  a  semi-infinite  elastic  medium.  The  results  obtained  are  for  both  coupled  and  uncoupled  cases 
and  are  compared  to  existing  analytical  solutions  for  their  accuracy.  A  comparison  is  also  given 
between  the  present  and  other  numerical  solutions. 


(i  A  KI:K  AMII  >  \S 


I.  ANALYTICAL  FORMl  LATION 


1.  Basic  Equations  and  Definitions 

Consider  a  medium  subjected  to  external  heating  and  assume  a  cartesian  coordinate  system 
(x,,i  =1,23).  It  is  assumed  that  the  medium  is  initially  at  uniform  temperature  T„,  which  will  be 
referred  to  as  the  reference  temperature,  and  the  state  at  this  temperature  as  the  reference  state.  It  is 
also  assumed  that  initially  the  medium  is  free  from  mechanical  deformations. 

The  instantaneous  absolute  temperature  is  denoted  by  T  and  the  difference  T  -  T„  defines  the 
instantaneous  iclative  temperature  Ad,  which  is  a  function  of  the  coordinate  x ,,  and  the  time  /.  Let 

T-  T0 


9 


(1) 


be  defined  as  the  temperature  change  per  unit  temperature  T0,  or  the  instantaneous  relative  tempera¬ 
ture  per  unit  temperature,  in  the  following  it  will  be  referred  to  as  the  temperature  0. 

In  the  formulation  of  this  study,  when  repeated  indices  in  equations  are  present,  the  summation 
convention  is  assumed.  We  now  define  a  vector  field  as  the  heat  displacement  vector  such  that 


9-**  - H 

9  dx,  "" 


(2) 


In  the  above  definition  the  temperature  9  represents  a  thermal  strain  similar  to  mechanical  strain. 
Note  that  H,  has  the  dimensions  of  displacement.  Thus,  there  is  a  one-to-one  correspondence  between 
heat  displacement-mechanical  displacement  and  temperature-strain. 

The  increment  of  entropy  s,  per  unit  volume,  for  the  case  of  zero  strain  is  given  by 

T-  T„ 

s-c-rr 

and  due  to  (1)  the  entropy  is  s  -  cfl,  where  c  is  the  heat  capacity  per  unit  volume  and  related  to  the 
heat  capacity  per  unit  mass  cf,  by  c  ~  pc„  with  p  as  the  mass  density  per  unit  volume. 
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By  introducing  a  quantity  a  as  the  heat  stress,  the  entropy  s  is  related  to  this  stress  by 

a  -  cT,fi  -  T0s  (3) 

Equation  (3)  can  be  considered  as  a  constitutive  relation  in  the  same  fashion  as  the  stress-strain  rela¬ 
tion. 


The  temperature  field  9  satisfies  the  thermal  diffusion  equation 


_3 

3x, 


kjj  (x, ) 

I 


B0(x„t) 

Bxj 


dO(X).t) 

Bt 


(4) 


wher  k,j  is  the  thermal  diffusion  coefficient  with  the  following  symmetric  property: 


V 


Using  Eq.  (2),  Eq.  (4)  is  expressed  as  follows 

3 

Bx, 
or 


h  i£.)  - 

3 

BH, 

3x, 

Bt 

BO 

'Bx, 


BH ) 
Bt 


and,  after  substituting  Eq.  (3)  in  the  last  equation,  yields 


B<r 

Bxj 


kn  ZZ--c'Ta 


BH, 

Bt 


or 


(5) 


Bxj  c  T°K° 


Bt 


Using  the  same  definitions  for  the  temperature  9  as  were  described  above,  and  under  the  assump¬ 
tion  of  small  mechanical  displacements  U,  for  every  particle,  the  kinematic  relations  are  given  by 


e,)  -  |  (U,j  +  Ujj)  -  /»*,  T0  l Cukl]-'9, 


0  “  Hu  —  —  Bn  eij’ 


(6) 

(7) 


5 


G  A  KERAMIDAS 


where  e,*  are  the  components  of  the  infinitesimal  mechanical  strain,  which  is  equal  to  the  difference 
between  the  total  strain  and  the  thermal  strain.  Similarly,  the  thermal  strain  6  is  equal  to  the  difference 
between  the  total  heat  strain  and  the  strain  induced  by  the  mechanical  deformation.  The  coefficients 
Cjjki  and  /)„  are  the  mechanical  and  thermal  moduli,  respectively,  which  are  symmetric  in  the  sense  that 

(-tiki  “  Cklij  =  C/ikl  ”  (-1, Ik  • 

@11  “  fi  If 

For  simple  one-dimensional  displacements  in  the  x-direction,  Eqs.  (6)  and  (7)  can  be  written  as 
follows 


e,  - 


BU, 

Bx 

Bx 


-  'fie. 


(8) 


The  constitutive  relations  for  coupled  thermal  deformations  can  be  expressed  in  terms  of  Eqs.  (6)  and 
(7)  as 

<r//  “  (-tiki  eki  ”  (-tiki  eu  ~  fiii  T0  9,  (9) 

a  -  cTffi  =  cTaH,  ,  -  /3„ToC,„  (10) 

where  <r„  are  the  components  of  the  mechanical  stress  tensor  and  v  is  the  heat  stress.  The  above  rela¬ 
tions  clearly  demonstrate  the  similarity  between  the  stresses  <rti  and  a. 


The  terms  involving  in  Eqs.  (6)  to  (10)  represent  thermomechanical  coupling  and  for  the 
uncoupled  case  one  usually  eliminates  the  coupling  effect  in  Eq.  (7)  but  retains  the  thermal  expansion 
effect  in  Eq.  (6),  since  the  temperature  variations  due  to  mechanical  deformations  are  negligibly  small 
in  many  cases.  The  temperature  then  is  evaluated  independently  of  the  strain.  The  momentum  equa¬ 
tion  is  given  by 


and  the  thermal  diffusion  equation  by 


PU, 


fi*; 

fix,  ' 


c2r0  A„/ir, 


8? 

fix, 


(11) 


(12) 
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where  p  is  the  mass  density  per  unit  volume,  A„  the  thermal  resistivity,  and  c  the  heat  capacitance  per 
unit  volume.  Body  forces  and  heat  sources  have  been  neglected  from  the  above  equations. 

An  important  feature  of  Eqs.  (11)  and  (12)  is  that  the  thermomechanical  coupling  term  does  not 
appear  explicitly  in  these  equations,  but  only  in  the  stress-strain  relations,  which  is  more  appropriate 
since  the  stress-strain  relations  are  those  which  describe  the  material  behavior.  The  equations  of  equili¬ 
brium  together  with  the  kinematic  relations  and  stress-strain  relations  will  be  used  next  to  derive  a  vari¬ 
ational  formulation  of  these  equations. 

2.  Variational  Formulation 

Choosing  as  unknown  variables  the  mechanical  displacement  U,  and  the  heal  displacement  H,  to 
describe  a  physical  system,  a  fundamental  variational  principle  is  derived  by  applying  the  principle  of 
virtual  work  to  the  equilibrium  equations. 

Consider  a  variation  8 U,  of  the  mechanical  displacement  U,  with  the  corresponding  variation  he,, 
given  by  Eq.  (6)  and  a  variation  hH,  of  the  heat  displacement  with  the  corresponding  variation  60  given 
by  Eq.  (7).  These  variations  are  assumed  to  be  virtual  displacements,  and  by  multiplying  Eq.  (11)  by 
hU,  and  Eq.  (12)  by  8 H,  ,  and  integrating  over  a  volume  v  of  the  medium,  one  obtains 

X  hU,dv  +  X„  (c2r0\„W,  -  v„)  hH,dv  *=  0  (13) 

Employing  the  divergence  theorem,  the  first  term  in  Eq.  (13)  becomes 

X,  (p&i  —  l/jdv  •»  JpUjSUidv  —  X,  (cr,j  SU,),jdv  +  J  a-,J(hU,),l  dv 

=  Jp  U,S  U,dv  -  X  (TIJnJ8U,dB  +  f  a-^e^dv  (14) 

=  X  P  Vfi  U,dv  +  X  Vifie„dv 

-  /, 

where  F,  is  the  surface  traction  applied  on  the  boundary  surface  B  of  volume  v. 
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Similarly,  the  second  term  of  Eq.  (13)  yields 

/  (c: 0Hj-  a. ,)6H,dv  -  X  ciT<>k,JHi  SH.J  v  -  f  (aS H,),,dv 

+  f  H,).,dv  (15) 

“  J(  c2T^ijHjh  Hjdv  +  cTffibOdv 

+ f  T^^&e.idv  -  fg  <r  n ,6 H,dB 

Hence,  Eq.  (13)  becomes 

X  p  Ufi  U,d\>  4-  X  c2T0kuHlbH,dv  +  £  cT^e  dv  +  X  o'  ifieijdv 

+  X 

=  X8  F,bU,dB  +  JganfiHtdB 
or 

X  pU,bU,dv  +  X  c^TokijHjBHjdv  +  X  Qw  ek,  6e„  dv 

+  X  cTo989dv  06) 

-  X,  [/-St/, 

Now,  if  one  considers  an  energy  function  which  is  expressed  as  the  product  of  stress  and  strain 

w  -  y  o-,,^,  +  y  (17) 

then  the  function  satisifies 

Bw  .  Bw  - 

— —  —  cr ji  and  — =-  —  cr  (18) 

8e„  1  B« 

where  9  is  the  total  heat  strain  and  equals 
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Alternatively,  writing  w  in  terms  of  strains 

"  -  y  (CUitekl  -  P0  T0e)e0  +  j  cTo0 
~  y  Cljkl  e,j  ek/  +  j  cTo07 


0  +  ~  PijCjj 


(19) 


8h>  -  Cijk,ellheki  +  cTo0h0 


(20) 


(21) 


A  potential  function  is  now  defined  as 

V-  y  /v  [Cilkle„ekl  +  cTo0}\  dv, 

which  for  isothermal  deformation  (0  -  0)  has  the  physical  meaning  of  a  strain  energy  function  and  for 
zero  strain  it  reduces  to  a  thermal  potential. 


(22) 


Variations  of  V  give 

8K-  /  Q,e(/8ewrfv  +  f  cT^hOdv. 

Hence,  Eq.  (16)  takes  the  form 

Jy  pU,BU,  dv  +  J  c'T'fiijHibH.dv  +8V  -  J*a  {f,8(/,  +  o6H,  n,]  rffi. 

The  second  term  of  Eq  (23)  is  related  to  dissipation  of  heat  energy  and  can  be  expressed  as  a  vari¬ 
ational  invariant  8Z>,  given  by 

8D-  c'TokAsH.dv  (24) 

and  Eq.  (23)  is  written  as 

J  pUt6Uidv  +  8 D  +  8K- fg  \FfiU,  +vn, dB  -  Q  (25) 

Equation  (25)  may  be  viewed  as  the  variational  form  of  the  coupled  thermal-momentum  diffusion 
equations,  with  Q  regarded  as  a  generalized  force. 
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3.  Generalized  Coordinates 


Consider  the  mechanical  and  heat  displacements  represented  by  the  following  forms 
U,  =  U^qx.q-i.  ....  q„,x„t) 

and  /=  1,2.3  (26) 

H,  =  H,(q,,q  2 . q„,  x„l) 

which  are  given  functions  of  the  space  coordinate  x,,  of  the  time  t  and  of  the  parameters  q,  .  These 
parameters  are  unknown  functions  of  time  and  will  be  considered  as  the  generalized  coordinates.  The 
variations  of  the  fields  U,  and  W,  with  respect  to  the  generalized  coordinates  q,  are  given  by 


...  _  at/,  -  3W,  /  -  1,2,3 

t.V,  -  —H.  SH,  -  —6," 


(27) 


The  variation  of  the  potential  V  is 


8  ^  ~ -  hq,  /  =  l,n 

dq, 


and  the  generalized  forces  Q,  have  the  form 


(28) 


/. 


»u,  OH, 

F' 


6 q,dB  -  Qfiq, 


and 


o-L 


p  Bo,  _  aw, 

F,  ~z —  +  an, 
dq,  '  dq, 


dB 


(29) 


Using  the  same  procedure  as  in  Lagrangian  mechanics,  the  first  term  in  Eq.  (25)  can  be  written  in 
terms  of  the  kinetic  energy  function  as  follows 

l  J>t4  Jj-  *q,dv. 

From  differential  calculus  one  derives 


dq,  dt 


_  Ut±  >JL 

dq,  dt  dq, 


(30) 
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and  from  Eq.  (26)  one  derives 


3  U,  bU, 

v‘  ”  h"  *' +  nr 


Taking  the  partial  derivative  of  the  above  with  respect  to  q , ,  yields 


a  U,  dU, 
a  <7,  3  qj 


and  with  respect  to  qk 

at/,  a2t/,  .  a;t/,  _  d  | at/, 

dqk  bqjbqk  +  bqkbi  dl  ( bqk 
Substituting  Eqs  (31)  and  (32)  into  Eq.  (30)  yields 


U  0  3tA  atA 

'  Bqj  d,  U'  bq,  '  3 q, 


If  the  kinetic  energy  is  expressed  as 


then  combining  Eqs.  (30)  -  (34)  one  derives 


f  pUibUfdv 

•r  v 


_d 

bK 

_  bt£ 

di 

bq, 

bq, 

(31) 


(32) 


(33) 


(34) 


(35) 


To  simplify  8  A  consider 


Hi 


bH,  .  aw, 
a <7,  *'  + 


and 


Thus 


bH,  9// 

a*,  "  3<7,  ’ 


8  Z> 


at? 

3$* 


«<7* 


II 
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where 


with  D  given  by 


BD 

Bqk 


f  c'T^H,  dVi 


B<lk 


(36) 


(37) 


Substituting  Eqs.  (28),  (29),  (35)  and  (36)  in  Eq.  (25)  one  obtains 


<A 

1 1>9, 

BK  |  BD  f  BV 
Bq,  Bq,  6  <7, 


0,69, 


(38) 


The  system  of  Eqs.  (38)  represents  the  Lagrangian  representation  of  the  governing  equations.  They 
constitute  a  system  of  n  differential  equations  corresponding  to  the  generalized  coordinates 


</,(/).  i  -  1,  n. 


Consider  now  a  special  case  where  the  displacements  are  approximated  by  a  linear  combination  of 
the  generalized  coordinates 

U,(Xj,i)  -  9*0)  /*,<*,)  U  “  >.2,3  (38) 

H((Xj,t )  =  Pk  0)  A,Or,)  *  *=  l.w 

(39) 

When  Eqs.  (39)  are  differentiated  with  respect  to  time  and  space  variables,  they  yield 

0,  -  qJk, ; 

H,  -  PkSkn  (40) 

O/, ,  “  9*  A /, 

Hij- pkgk,.j.  <4I) 

Then  the  strains  e„  and  9  are  given  by 

e„  -y  9t  (/*,.,  +  fklJ 
«  -  />*#*,.,  ~  ~~  e' 

12 
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The  potential  energy  Kean  be  expressed  as 


V 


1  I  2  1  3 

J  viMi  ~  vfjQiPj  +  y  v^pj, 


(43) 


where 


v« 


"X 


To 


| Qjkl  +  /8<A,|  (fklJ  +  W 

vk!  -  X  +  /^,)  ginmdv, 

vki  *  cT'o  X,  &,./  gijj- 

Similarly,  the  dissipation  and  kinetic  energy  functions  are 


D  **  y  dkiPkPi- 


where 


K'y  mkiQk9i- 

dki  -  X  c'ro  Kugk,gi,clv, 
”<kl  «  X  Pfkifhdv- 


The  corresponding  Lagrangian  equations  are  then  given  by 

m„iij  +  v,^  -  v2p,  -  FF, 

d„Pj  +  v,3,p,  -  v2"9/  -  FQ, 


4y  »  1  ,n 


(44) 


(44) 


(45) 


(46) 


(47) 


where  FF,  and  FQ,  are  the  mechanical  and  thermal  forces  corresponding  to  the  generalized  coordinates 
q,  and  p,  respectively,  and  given  by 

FF'  "  X  a i^kfudB.  FQ,  -  fB  vrt'g.jdB.  (48) 

Eqs.  (47)  represent  a  system  of  2  n  differential  equations  for  the  unknown  field  parameters  q,  and  p, 
which  describe  the  mechanical  and  heat  displacement  fields.  The  solution  of  this  system  then  gives  the 
thermal  and  mechanical  deformations. 

Note  that  the  thermomechanical  coupling  is  included  in  the  coefficient  v2.  For  the  uncoupled 
case,  the  term  v,2,g,  is  neglected  and  the  second  equation  may  be  solved  independently  without  the 
knowledge  of  q,. 
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II.  FINITE  ELEMENT  FORMULATION 

The  formulation  in  this  study  includes  two  kinds  of  one  dimensional  finite  element  models:  a)  a 
linear  element  (LE),  with  the  minimum  number  of  degrees  of  freedom,  and  b)  a  higher-order  element 
with  four  degrees  of  freedom,  two  for  each  nodal  point,  which  will  be  referred  to  as  the  cubic  element 
(CE).  The  disadvantage  of  the  linear  element  is  that  the  value  of  strain  throughout  the  element 
depends  only  on  time.  To  achieve  good  results,  one  has  to  use  a  large  number  of  elements,  particularly 
for  problems  which  involve  transient  behavior  and  discontinuities. 

It  should  be  pointed  out  that  existing  finite  element  analyses  are  usually  based  on  the  approxima¬ 
tion  of  displacements.  For  pure  thermal  diffusion  problems  the  temperature  may  be  treated  as  a  gen¬ 
eralized  displacement  to  apply  finite  element  approximations.  However,  when  other  types  of  behavior 
are  involved,  such  as  mechanical  deformations,  similar  treatments  fail  to  yield  a  compatible  formula¬ 
tion.  This  is  the  main  reason  that  the  present  formulation  is  given  in  terms  of  the  heat  displacement 
field. 

1.  The  Linear  Element  Analysis 

The  displacement  fields  U,  and  H,  are  approximated  by 

U(x,i)  -  o i (/)  +  ai(t)x 

and 

H(x,t)  -  *,(/)  +  Mf)x  (49) 

The  conditions  at  the  two  nodal  points  are 

at  x  -  0:  U  -  Uu  H  -  H\ 
at  x-  I.  U-  U2,  W  -  Hj 
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where  C £/(,  C/2)  and  (H\,H2)  are  the  nodal  displacements.  Applying  the  above  conditions  in  Eq.  (49) 
one  obtains 

£/-(/,  +  (U2-  Ut)j 

//-//,  +  (Hj-HOj 

or,  in  an  alternate  form 


U 


U,  +  jU2 


H 


(50) 


One  can  identify  the  shape  functions  /„,  g0  and  the  generalized  coordinates  q,  and  p,  as 


“  /it  "  (1  -  ■yj.  fu-/i2“  -y; 

Q\m‘  U\,  q2~  £/2; 

Pi  -  p2  -  H2. 


(51) 


The  only  non-vanishing  strain  component  is 


(52) 


and  the  stress  components  are  then  expressed  as 

a a  “  C/yii^n  ~  PijT 

If  an  isotropic  elastic  material  is  considered,  shear  stresses  are  zero  and  the  axial  stress  com¬ 
ponents  are 

& u  “  (X  +  2ft)eii  - 

<5,) 

(r22  ~  *33  ”  “  ®|’ 
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The  kinematic  relation  gives 


-y  (H2~  «,)  -  W2~  (/,)|. 


The  matrix  coefficients  from  Eqs.  (44)  and  (45)  can  be  evaluated  and  the  results  are 


1  V  .  ,  t*  '  o 

v,1,  -  x  f  2 m  + -  a„; 

c 


V*  -  PT0a,r 
v,;  -  cT0air 
cT0 

d"  ~  ~1T  b"' 

m„  -  pb„. 


I  1  -1  A 

where  a„  -  I,  ,  —  and  b„ 
displacements  are 


2  1  Al 

I  j  Then  the  governing  equations  for  calculating  the  nodal 


pb„U,  +  A  +  2p  4-  a„U,  -p  T0a„H,  -  FF„ 


—  b.jH,  +  cT0a,jH,  —  (3T0atl(Jj  =  FQ„ 


FF,  -  A(3  Ttfr ,n 


FQ,  -  AcT<fl,n. 


where  <r,  is  equal  to  <rn  for  the  nodal  point  1.  For  the  particular  case  of  uncoupled  thermoelasticity  , 
the  temperature  variation  due  to  coupling  is  very  small  and  the  /3r0o„l/,  term  in  the  second  equation 
may  be  neglected. 


16 


NRL  MEMORANDUM  REPORT  4310 


2.  The  Cubic  Element  Analysis 


Assume  the  following  approximations  for  the  displacement  fields 


U  -  <j0  +  a\x  +  aix 2  + 


H  -  +  b\x  +  b2x 2  +  b}x}\ 


e  “  TT  "  +  2o2x  +  3a}x2 


—  — -fie  —  b]  +  2/)2x  +  ibyx2  — -fie 

fix  c  '  2  c 


The  conditions  at  the  two  nodal  points  are 


at  x  -  0:  V  -  U\,  e  -  et.  H  -  Hx  and  f»  -  0,  -  —  fie , 
at  x  -/:£/—  £/2,  e  —  e2,  H  —  //2  and  9  02  -  -jfie2 

where  (f/,,f/2)  and  (HitH2)  are  the  nodal  displacements  and  (C|,e2),  (0,,02)  are  the  nodal  mechanical 
and  thermal  strains.  The  coefficients  a,  and  b ,  are  identified  as 

oo~  U\>  °\m‘  e\ • 

a2~-  j  l(e2  +  2 e,)l  -  3 (U2-  t/,)]. 

«j-  4l(e2  +  c,)/-  2(t/2-t/,)l.  (60) 

bo“  H\,  b\  “  8), 

b2  -  -  -y  t(0i  +  20,)/  -  3(//2  —  //,)). 

*3-  "jj  l(02  +  0i)/  —  2(H2—  H\)\. 


Hence,  the  displacement  and  strain  components  have  the  form 

U  m  f\\U\  +  /l2*l  +  f li^l  +  /|4*2- 

e  -  011^1  +  *i2*t  +  g\iUi  +  g\*ei. 

H  "  f\\H\  +  /i^i  +  ftjH2  +  /m02. 

®  “  g i|H|  +  /i^i  +  guW2  +  /M02  -  —  fie. 
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One  may  identify  the  shape  functions  /1(  as 


/i.= 


2x3 
/J  • 


/ 12 


X 


/l4 


and  gu  is  given  by 


(62) 


Si/ 


9/i/ 

9x 


The  generalized  coordinates  q,  and  p,  can  be  identified  as 


<7 1  -  */i.  <72  =  *1.  <73  =  i/2.  <74=  <”2. 
Z7!  “  //|,  =  <Z|,  =  //j,  Z>4  =  #2- 


The  expression  for  v*,  </„,  and  can  be  written  as 


(A  +  2/z )  + 


P2To 


&ii’  v ii  PT0a„,  v,',  *  cTqOu, 


cTn 


b„.  m„  ~  pb„, 


(64) 


where  0,,  and  6„  are  given  by 


and 


A 

30/ 


_Al_ 

420 


36 

3/ 

3/ 

4/2 

-36 

-3/ 

3/ 

-/2 

156 

22/ 

22/ 

4/2 

54 

13/ 

-13/ 

-3/2 

-36  3/ 

-3/  -/2 
36  -3/ 
-3/  4/2 


54  -13/ 
13/  — 3/2 
156  -22/ 
-22/  4/2 


(65) 


(66) 
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The  governing  equations  in  terms  of  the  generalized  coordinates  are 


pb„Q,  + 


(A  +  2m)  + 


P2T0 


a„Q>  ~  PToa.'P,  =  FF, 


-p-  Ki><  +  cToa„p,  ~  0  top,^,  -  FQ, 

with  the  generalized  forces  FF,  and  FQ,  given  by 

FF)  ”  Af3T(fr,n,  i  -  1,3, 

FQ,  —  AcT<fl,n,  /  —  1,3, 
and 


(67) 


(68) 


/T2  -  -  F02  -  f(?4  -  0. 

The  differential  Eqs.  (67)  constitutes  a  system  of  eight  equations  for  the  nodal  displacements  and  nodal 
strains. 


3.  The  Overall  Problem  and  Boundary  Conditions 

In  the  previous  sections,  the  equations  for  the  basic  linear  and  cubic  element  models  were 
derived.  Consider  now  a  one-dimensional  case  with  the  medium  modeled  by  either  of  the  element 
models.  Since  the  heat  and  mechanical  displacements  are  compatible,  the  connection  of  the  elements 
follows  the  rules  and  the  assembly  process  analogous  to  the  one  used  in  the  matrix  analysis  of  struc¬ 
tures.  The  most  commonly  used  assembling  method  is  the  direct  stiffness  method,  according  to  which 
the  overall  stiffness  matrix  results  from  the  stiffness  matrices  of  the  individual  elements  by  simple  addi¬ 
tion  at  the  nodal  points. 

For  the  elements  introduced  in  this  analysis,  the  matrix  equation  for  the  element  can  be  put  in 
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Pb.j  0  q,  + 

0  oj  />,} 


0 

0 


0 

cT0 


Q, 

P, 


+ 


(X  4-  2 fi )  + 


/s2r0 


~P  Toa„ 


(3T0  o,j 
cT0  a„ 


(69) 


Then  the  assembled  matrix  equation  has  the  following  schematic  form 


Mi)  o 

( a) 

t 

0  0 

!<?}] 

i 

IA}) 

Ia4) 

(<?)  _ 

0  0 

|i/>) 

0  [Aj] 

\p)\ 

“T 

[AS] 

lAb) 

{Ip! 

{FQ) 

where  the  square  matrices  [A,]  are  derived  through  the  direct  stiffness  method.  The  vectors  Ig)  and 
[p)  represent  nodal  degrees  of  freedom  and  the  force  vectors  | FF)  and  (FQ!  have  zero  components  at 
the  connecting  nodal  points  and  are  non-zero  on  the  boundary.  The  solution  of  the  above  matrix  equa¬ 
tion  gives  the  nodal  displacements.  Then  the  distribution  of  stresses  can  be  evaluated  by  the  stress- 
displacement  relations  and  the  kinematic  relations. 


The  formulation  of  the  overall  problem  is  not  complete  unless  boundary  conditions  are  taken  into 
consideration.  The  procedure  is  identical  to  the  one  in  traditional  matrix  structural  analysis,  where  par¬ 
titioning  of  the  equations  eliminates  the  singularity  of  the  matrix  coefficients. 
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III.  APPLICATION  OF  THE  FINITE  ELEMENT  FORMULATION 
TO  THE  COUPLED  DIFFUSION  PROBLEM 

1.  Problem  Formulation 

In  applying  the  derived  finite-element  formulation,  the  problem  of  the  linear  elastic  half-space 
subjected  to  a  time  dependent  temperature  change  on  its  traction-free  boundary  plane,  is  considered. 
The  initial  boundary  value  problem  assuming  a  sudden  heating  of  the  boundary  plane  is  known  as 
Danilovskaya’s  problem  and  an  extension  of  this  problem  is  given  by  Sternberg  and  Chakravorty  in 
Ref.  [20]. 

The  specific  notations  and  descriptions  of  the  problem  are  depicted  in  Fig.  1.  Let  ( x,y,z )  be  the 
Cartesian  coordinate  system  and  consider  an  elastic  medium  occupying  the  half-space  x  >  0,  with  the 
boundary  plane  at  x  =  0  assumed  to  be  free  of  traction  at  all  times.  For  uniform  boundary  heating,  it 
can  be  assumed  that  the  medium  is  restricted  to  a  uniaxial  motion  and  hence 

£4-  Ux(x.i). 

(71) 

V,  =  U.  =  0. 

For  the  displacements  in  Eq.  (71),  the  Cartesian  components  of  shear  stress  vanish  identically  and  the 
stress-strain  relations  for  a  linear  material  are  given  by 

cr  r  =  (X  +  2fi )  ex  —  /3  ( T  —  Tq), 

(72) 

<tv  -  <r2  -  kex  -  p(T  -  T0). 

where  (<r,,<rv,<r.)  are  the  components  of  normal  stress,  X,  n  are  the  material  constants,  p  is  the  ther¬ 
mal  modulus  and  T(x.i)  is  the  time  dependent  temperature  field  and  T0  -  T(x,  0). 

If  the  body  forces  are  absent,  then  the  momentum  equation  reduces  to 


Ci.  A  KtRAMlOAS 


in  which  p  denotes  the  constant  mass  density,  and  the  temperature  field  is  governed  by  the  thermal 
diffusion  equation 


BT  _  B2T 
dt  ~  *  dx 2 

where  «  =  k/c  is  the  thermal  diffusivity.  The  initial  condition  for  the  temperature  is 


(74) 


T(x,0)  =  r0. 

while  the  uniform  ramp  heating  yields  the  boundary  conditions 


(75) 


T  —  T0  —  0  for  t  <  0; 

T  ~  r0=  —  (T,  -  T0)  for  0  ^  f  <  r0;  (76) 

1 o 

=*  T{  -  T0  for  t  >  r0; 

where  r0  represents  the  time  required  for  the  temperature  to  reach  a  constant  value  T\. 


If  the  medium  is  initially  at  rest,  then  the  initial  conditions  for  the  displacement  are 


bUx 

UM  0)«0,  -^U,0)-0, 

dt 


while  the  stress  boundary  condition  is 


(77) 


o-x(0,r)  =  0,  t  >  0. 


(78) 


To  satisfy  the  regularity  requirements,  one  also  has 


(T(x,t)  -  T0),  Ux{x,t),  <rr(x,t),  tr  ,,(x,r)  —  0,  as  x  — 


(79) 


At  this  stage  it  is  expedient  to  relate  the  dimensionless  variables  to  the  physical  ones  as  follows: 


f  -  —  T  - 

K  K 


O’ I  - 

*2‘ 


P(T,  -  T0) 

- - Zi - 

/3(T,  -  T0) 


,  U  - 


H  - 


a  (A  +  2 p )  £/x 
W-  t0)k  ' 
o  7-p//, 
k(T,-  T0)' 


(80) 


9  T~  T0  ^  T0 


#,  r,-r0  r,-r0' 

where  k  *=  A/c,  -  (A  +  2pt)/p  and  /3  a(3A  +  2p),  with  a  being  the  thermal  coefficient  of  linear 


expansion.  H  is  the  dimensionless  heat  displacement  in  the  present  formulation. 
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Eqs.  (73)  and  (74)  can  be  written  as 


and 


df  “  8r2' 


6£  =  3^0 
3r  3f 2  ‘ 


The  initial  conditions  are 


and  the  boundary  conditions  are 


U(f.  o)« 

OT 


Type  1.  Sudden  heating  of  the  boundary 


cri  (0,  r)  “  0  for  all  r; 

|0  r  <  0 

0(°,t)«  j  T  ^  0; 


Type  2.  Ramp-heating  of  the  boundary 


o-|(0,  r)  =  0  for  all  r; 


0(0,  t)  = 


0 


T  <  0 


t/tq  0  <  t  <  t0. 

I  r0  <  T 


(81) 


(82) 


(83) 


(84) 


Applying  the  finite-element  formulation  developed  in  the  previous  chapters,  the  equations  for  the 
two  element  models  can  be  written  in  the  following  dimensionless  forms  given  below 

Linear  element  (LE)  model: 


(6,1(01 
(0)  (01 


lift 

lift 


(01  (01 

(01  (6,1 


u» 


+  6h>2 


(1+8)  [a, 1  ~(o,l 
-8  [a,]  [a,] 


n 

—  6w 

IFF) 

\lH)\ 

IFQ) 

(85) 
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Here,  w-x/al  with  /  as  the  element  length  and  FF,  -  (FF,)l[Ac(Tx  -  r0)l, 
FQ,  -  (FQ,)/[Ac(T\  -  r0)l.  The  two  constant  matrices  [a,]  and  [8,1  are  given  by 


la,] 


1  -1 
-1  1 


and 


Cubic  element  (CE)  model: 


(86) 


lb,  ]  [0] 

[0]  [0] 


If/) 

Ve) 

|W> 

l«) 


[0]  [0] 
f0]  (dj 


10) 
le) 
I H) 


+  14tv2 


(1+8)  [ac]  -  [a,  ] 
-8  [a(]  [a,.] 


M 

{£/) 

\e) 

«=  420  tv 

10} 

1 H) 

1  FQ) 

\9\ 

(0) 

(87) 


Here  the  parameter  8  is  defined  as  the  thermomechanical  coupling  parameter,  given  by 


8 


and  the  matrices  [af]  and  [6f]  are  given  by 


PT0 

c(A  +  2fi)  ' 


(88) 


36 

-36 

3/  tv 

3/  tv 

-36 

36 

-3/t v 

-3/ tv 

3/ tv 

—3/ tv 

4/  tv2 

—  1/  tv2 

3/ tv 

-3/  tv 

—  1/  tv3 

4  tv2 

lbc) 


156  54  22/ h>  -13/ tv 

54  156  13/ tv  -22/ w 

22/ w  13/ tv  4/  tv2  -3/ tv2 

-13/ tv  -22/ w  —3/  tv2  4/ tv2 


(89) 


The  assembly  of  the  above  equations  for  the  overall  problem  and  their  modifications  due  to  the  boun¬ 
dary  conditions  was  coded  in  the  computer  program. 
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In  the  case  of  the  (LE)  model,  after  solving  for  the  displacements,  the  temperature  for  the  /  th 
element  can  be  obtained  from  the  relation 

0,  -  -  H,\  -bw\Ut+x-  U,].  (90) 

and  the  stresses  for  the  i  th  nodal  point, 

o-i,  =  w(Ui+l  -  Ut)  -  0,.  (91) 

cr  2,  =  v<r\j  —  (I  —  2y)0(. 

For  the  (CE)  model,  the  solution  of  the  equations  will  directly  give  the  nodal  displacements  and 
strains,  and  thus  the  temperature  and  strain  distributions  are  obtained  directly  from  the  solution  of  the 
matrix  equation.  Eqs.  (85)  through  (91)  can  be  applied  to  the  coupled  and  uncoupled  cases  by  simply 
setting  8  ^  0  or  8  =  0,  respectively.  These  equations  can  be  used  for  solving  the  problem  stated  pre¬ 
viously  in  this  section. 

2.  Numerical  Results 

A  numerical  solution  of  Eqs.  (85)  and  (87)  for  the  finite-element  analyses  is  obtained  by  using  a 
third  order,  backward  finite-difference  scheme.  A  closed  form  analytical  solution  is  given  in  Ref.  [20] 
for  the  uncoupled  equations  and  it  will  be  used  here  to  represent  the  exact  solution.  An  analytical  solu¬ 
tion  of  the  coupled  case  is  given  by  Nickell  and  Sackman  in  Ref.  [16],  for  the  ramp-heating  type  of  the 
boundary  condition. 

For  the  finite-element  solution  the  total  number  of  elements  used,  was  TNE  “  30  for  the  (LE) 
model  and  TNE  =  20  for  the  (CE)  one.  The  temperature  boundary  condition  is  characterized  by  the 
parameter  t0  as  follows: 

Type  1 .  Sudden  heating,  r0  —  0 

«<o,r>  -  {»  ;  <  »  <v.M> 
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Type  2.  Ramp-heating,  t0  -  1.0 

0  r  <  0 

0(0, r)  -  t/t0  0  <  t  <  1.0  (V.2.2) 

1  r  >  1.0 

Two  values  of  the  coupling  parameter  were  used,  8  «  0  for  the  uncoupled  case  and  8  —  1  for  the 
coupled  one.  The  time  step  sizes  used  for  solving  the  matrix  differential  equations  were  Dt  —  0.005 
for  the  (LE)  and  Dr  —  0.01  for  the  (CE).  For  evaluating  the  lateral  stress  <r2,  the  value  of  v  —  1/3 
was  used  for  Poisson's  ratio. 

The  results  obtained  are  illustrated  in  Figs.  (1)  through  (12)  for  both  element  models  and  they 
are  compared  to  the  exact  solution.  For  the  case  of  r0  -  0,  8  —  1,  numerical  data  for  the  exact  solu¬ 
tion  were  not  available,  thus  only  the  results  from  the  present  formulation  are  given.  Figs.  1  to  4  dep¬ 
ict  the  temperature,  mechanical  displacement  and  stresses  as  a  function  of  time  at  {  —  1.0  and  for  r0  — 
0.  The  results  for  the  temperature.  Fig.  1,  show  a  very  good  agreement  with  the  exact  solution  for 
both  element  models.  For  the  displacement,  Fig.  2,  the  (CE)  model  gives  less  error  than  the  (LE), 
with  very  small  instability  around  the  point  of  discontinuity.  For  the  stresses,  Figs.  3  and  4,  again  the 
results  for  the  (CE)  are  better  than  those  of  the  (LE).  The  instability  around  the  discontinuity  point  r 
=  1.0,  is  due  partly  to  the  numerical  scheme  and  partly  to  the  small  number  of  elements  between  f  “ 
0  and  £  =  1.0.  This  instability  can  be  corrected  by  refining  the  mesh  around  the  point  of  discontinuity. 

The  results  for  r0  “  1.0,  are  illustrated  in  Figs.  5  to  8.  The  temperature  as  a  function  of  time  at 
£  -  1.0,  is  given  in  Fig.  5,  the  displacement  in  Fig.  6,  and  the  normal  and  lateral  stresses  in  Figs.  7  and 
8,  respectively.  As  one  can  see  from  these  figures,  the  results  obtained  through  the  CE  solution  show 
excellent  agreement  with  the  exact  solution  even  around  the  discontinuity  point.  For  the  (LE)  solu¬ 
tion,  the  results  for  the  displacement  and  stresses  show  an  error  which  is  due  to  the  numerical  approxi¬ 
mation  and  can  be  corrected  by  refining  the  mesh  around  the  point  f  ™  1.0.  In  order  to  show  the  pro¬ 
pagation  of  discontinuities  in  the  half-space,  the  temperature  is  plotted  in  Fig.  9,  the  displacement  in 
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0.0  0.4  0.8  1.2  1.6  T  -  2.0 

Figure  6  —  Temperilure  distribution  et  f  •  I.  is  ■  function  of  time,  with  t0  -  1  0.  dynimic  solution 
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Figure  8  —  Normal  stress  distribution  at  {  -  I.  as  a  function  of  time,  with  r„  -  |  0,  dynamic  solution 


30 


G.  A.  KERAM1DAS 


Fig.  10,  and  the  normal  stress  in  Fig.  11  as  a  function  of  the  space  coordinates,  for  t0  ”  0.25  and  t  — 
1.0,  r  —  2.0.  These  figures  show  how  the  thermal  induced  wave  propagates  into  half-space.  As  one 
can  see  from  these  figures,  the  (CE)  gives  excellent  agreement  with  the  exact  solution,  and  the  LE 
gives  satisfactory  results  but  with  a  small  instability  around  the  discontinuity  points. 

SUMMARY  AND  CONCLUSIONS 

A  unified  formulation  for  the  coupled  thermal-momentum  diffusion  equations  was  introduced  in 
this  report  and,  based  on  this  formulation,  two  finite  element  models  were  developed  for  the  purpose  of 
solving  problems  on  thermal  and  mechanical  deformations. 

The  introduction  of  a  new  quantity,  defined  as  heat  displacement,  is  the  basis  for  a  unified  presen¬ 
tation  of  the  governing  equations  with  one-to-one  correspondence  between  thermal  and  mechanical 
quantities.  This  presentation  is  used  to  develop  a  successful  displacement  formulation  for  the  finite  ele¬ 
ment  analysis.  The  resulting  matrix  differential  equations  for  the  finite  element  models  are  given  in 
terms  of  nodal  displacements  and/or  strains.  Due  to  the  nature  of  this  formulation,  boundary  condi¬ 
tions,  given  in  terms  of  temperature,  strain,  mechanical  displacement  and/or  heat  flow,  can  be  easily 
handled. 

In  the  present  formulation  a  thermal  force  was  also  introduced,  for  which  one  should  point  out  its 
significance  as  a  boundary  force.  The  regularity  condition  for  Danilovskaya’s  problem  requires  that 
0—  0  as  £—•«>.  Since  the  last  nodal  point  in  the  finite  element  approximation  of  the  half-space 
represents  infinity  one  should  impose  the  above  condition  at  this  point,  and  then  the  thermal  force  is 
zero  due  to  zero  temperature.  This  assumption  is  not  the  correct  one  since  the  temperature  at  the  last 
nodal  point  will  increase  as  the  thermal  wave  propagates  and  the  time  r  increases.  If  we  consider  the 
last  nodal  point  as  a  boundary  point  and  the  thermal  force  as  a  boundary  force  which  is  equal  to  the 
temperature  at  that  point,  then  the  conditions  of  the  boundary  point  are  properly  adjusted.  The  impor¬ 
tance  of  this  thermal  force  is  shown  in  Fig.  12  ,  where  the  temperature  is  given  as  a  function  of  (  for 
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1.0  2.0  (  *  3.0 

Figure  11  —  Displacement  distribution  as  a  funcuon  t.  dynamic  solution  with 
r o  *  0  25.  as  a  function  of  time,  with  r0  •  0  0,  dynamic  solution 
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Figure  12  —  Normal  stress  distribution  as  a  function  f .  dynamic  solution  with  r0  -  0.25 
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the  three  different  times  r  and  for  r0  —  0.0.  The  results  illustrated  show  how  the  finite  element  solu¬ 
tion  behaves  by  neglecting  or  retaining  the  thermal  forces.  When  the  thermal  force  is  neglected  the 
solution  is  lower  than  the  exact  solution  for  {  >3  and  as  r  increases,  but  this  is  not  true  if  the  thermal 
force  is  retained  in  the  formulation. 

A  solution  of  the  matrix  differential  equations  can  be  obtained  by  any  standard  integration  tech¬ 
nique.  Here  a  third  order  backward  finite-difference  scheme  was  adopted  for  the  solution. 

The  two  finite  elements  developed  in  the  present  study  were  used  successfully  to  model  the  pro¬ 
pagation  of  thermally  induced  waves  in  a  semi-infinite  medium,  with  prescribed  conditions  for  the 
stress  and  temperature  on  the  boundary  plane.  For  the  temperature,  the  boundary  conditions  are  of 
two  types;  one  is  sudden  heating  and  the  other  is  ramp  heating  is  the  boundary  plane.  Results  obtained 
by  the  present  finite  element  analysis  are  for  the  cases  both  of  coupled  and  uncoupled  behavior.  These 
are  compared  to  existing  analytical  solutions.  The  cubic  element  model  gives  an  excellent  representa¬ 
tion  of  the  exact  solution  and  the  linear  element  also  gives  satisfactory  results  within  acceptable  limits 
of  accuracy. 

A  comparison  of  the  present  solution  with  existing  numerical  solutions  of  similar  problems  is 
appropriate  in  order  to  evaluate  the  efficiency  of  the  present  formulation.  Two  of  the  numerical  solu¬ 
tions  which  solve  the  same  thermoelastic  problem  are  given  by  Nickel!  and  Sackman  in  [16]  and  by 
Oden  [17]. 

The  solution  given  by  Nickell  and  Sackman  is  based  on  the  Ritz  method  and  their  results  give 
satisfactory  agreement  with  the  exact  solution.  In  order  for  the  authors  to  achieve  a  stable  solution, 
with  small  error,  they  used  a  very  fine  mesh  for  the  characteristic  length  L,  (TNE  —  47,  NE  —  20), 
and  a  time  step  of  0.005-0.01.  The  computing  time  for  such  solutions  was  about  five  minutes.  The 
present  solution,  with  a  less  fine  mesh,  produces  results  of  similar  accuracy  and  with  less  computing 
time  required  to  solve  the  same  problem.  For  example,  the  linear  element  model  solution  with  TNE  *■ 
30  requires  about  70  seconds  and  the  cubic  one  with  TNE  -  20,  requires  about  100  seconds.  Thus, 
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the  present  solution  is  more  efficient  for  the  same  degree  of  accuracy.  Another  advantage  of  the 
present  formulation  over  the  one  given  by  Nickel!  and  Sackman  is  that  it  can  be  used  for  any  geometric 
configuration  and  with  any  kind  of  boundary  condition. 

The  second  solution,  given  by  Oden,  is  based  on  the  finite  element  method.  In  reference  117], 
Oden  gives  results  only  for  the  temperature  and  displacement.  For  some  of  his  results,  Oden  used  up  to 
TNE  **=  48  with  NE  =  10  in  order  to  achieve  a  satisfactory  solution.  Comparing  the  two  solutions  it  is 
apparent  that  the  one  proposed  here  requires  a  smaller  number  of  elements  for  the  same  accuracy  as 
Oden’s.  Oden  also  mentions  that  his  results  for  the  temperature  are  lower  than  the  exact  solution  for 
£  >  3  and  as  the  time  r  increases,  and  he  suggests  an  increase  of  the  characteristic  length  L  in  order  to 
obtain  improved  results.  This  kind  of  correction  becomes  unnecessary  to  this  formulation  since  this  is 
corrected  by  the  thermal  force  of  the  last  nodal  point  as  was  shown  in  Figure  12. 

In  conclusion,  the  unified  form  of  the  present  formulation  should  be  emphasized,  as  well  as  its 
efficiency  on  handling  various  types  of  boundary  conditions.  The  superiority  of  the  cubic  element 
model  over  the  linear  type  should  be  also  noted.  Both  element  models  can  be  used  for  solving 
diffusion  problems  and  the  choice  between  them  depends  on  the  particular  needs  of  a  problem.  An 
extension  of  this  formulation  to  other  problems  or  to  two  dimensions  is  a  next  step  in  applying  the 
finite  element  model. 
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